% % ADMM-based OFDM Radar Waveform Design
% % Objective: Minimize spectral ISL via unimodular sequence
% % Constraints: Time-domain oversampled near-unimodularity, and unimodular constraint

close all;clear

%% Parameters
N = 128;                % Sequence length
M = 4*N;                % Oversampling factor
rho = 5e-1;             % ADMM penalty parameter
maxIter = 1000;
epsilon = 1e0;          % Time-domain unimodularity tolerance

%% Initialization
u = exp(1j*2*pi*rand(N,1));           % Initial unimodular sequence
u_init = u;
F1 = dftmtx(2*N);                      % DFT for frequency domain (2N points)
F = F1(:,1:N);
F1 = dftmtx(M); 
F1 = (M / N) * F1' / M;
% % index4 = zeros(1,M-3*N);
% % index4(1:1:end) = 4*(1:N) - 2;
% % index4(2:2:end) = 4*(1:N) - 1;
% % index4(3:3:end) = 4*(1:N);
T = F1(:,1:N) * dftmtx(N);
% T = ifft(eye(M,N));                  % Oversampled time-domain operator (e.g. interpolation matrix)
z1 = F * u;            % Frequency domain sample (zero-padded)
z2 = T * u;                          % Oversampled time-domain samples
z3 = u;                              % Copy for unimodular projection
y1 = zeros(2*N,1);
y2 = zeros(M,1);
y3 = zeros(N,1);

osu = ifft(cat(1, fft(u), zeros(3*N, 1)));
acf = conv(conj(flip(osu)), osu);
acfmag = abs(acf);
acfdb = 20*log10(acfmag / max(acfmag));
figure(1);hold on;legend();
title('过采样时的自相关（dB）')
plot(acfdb, 'DisplayName', '起始自相关函数')

figure(2);hold on;legend();
title('过采样后的包络');
plot(abs(T * u), 'DisplayName', '起始过采样包络');

figure(3);hold on;legend();
title('恒模序列过采样后功率谱')
plot(abs(fft(T * u)), 'DisplayName', '优化前');

figure(4);hold on;legend();
title('恒模序列过程采样的ACF模值')
plot(acfmag, 'DisplayName', '优化前')

figure(5);
subplot(121)
plot(abs(z2))
ylim([0 2])
subplot(122)
plot(abs(z1))

%% ADMM Iteration
rho1 = 1*rho;

rho2 = 10000*rho;
rho3 = rho;
lamda1 = 6;

% A = rho1*(F'*F) + rho2*(T'*T) + rho3*eye(N);
A = rho1*(F'*F) + rho2*(T'*T);
desired_mag = sqrt(N);
  
for k = 1:maxIter
    % === u-update ===
    
    % b = rho1*F'* (z1 - y1/rho1) + rho2*T'*(z2 - y2/rho2) + rho3*(z3 - y3/rho3);
    
    % A = rho2*(T'*T) + rho3*eye(N);
    % b = rho2*T'*(z2 - y2/rho2) + rho3*(z3 - y3/rho3);

    % A = rho1*(F'*F) + rho3*eye(N);
    % b = rho1*F'* (z1 - y1/rho1) + rho3*(z3 - y3/rho3);

    b = rho1*F'* (z1 - y1/rho1) + rho2*T'*(z2 - y2/rho2);
    
    u = A \ b;

    % === z1-update (spectral magnitude matching) ===
    Fu = F * u;
       
    % psi = angle(Fu);
    % target = desired_mag * exp(1j * psi);
    % z1 = (rho1*Fu + rho1*y1 + 2*N*target) / (rho1 + 2*N);

    psi = angle(Fu + y1 / rho1);
    target = desired_mag * exp(1j * psi);
    % z1 = (rho1*Fu + rho1*y1 + 2*N*target) / (rho1 + 2*N);
    z1 = (2 * lamda1 * desired_mag + rho1 * abs(Fu + y1 / rho1)) / (2 * lamda1 + rho1);
    z1 = z1 .* exp(1j * psi);
    y1 = y1 + rho1 * (Fu - z1);

    % Optional: monitor convergence
    if mod(k,1) == 0
        fprintf("Iter %d: ISL cost = %.4f\n", k, sum((abs(Fu) - desired_mag).^2));
    end

    % === z2-update (project to near-unit circle with energy constraint) ===
    Tu = T * u;
    c = Tu + y2/rho2;
    theta = angle(c);
    z2 = exp(1j * theta);
    if norm(z2 - c)^2 > epsilon
        % Relaxation to satisfy energy constraint
        alpha = sqrt(epsilon) / norm(z2 - c);
        z2 = c + alpha * (z2 - c);
    end
    
    figure(5);
    subplot(121)
    plot(abs(Tu))
    ylim([0 2])
    subplot(122)
    plot(abs(Fu))

    osu = ifft(cat(1, fft(u), zeros(3*N, 1)));
    acf = conv(conj(flip(osu)), osu);
    acfmag = abs(acf);
    acfdb = 20*log10(acfmag / max(acfmag));
    figure(6);hold off
    plot(acfdb)
    ylim([-60 0]);
    grid on;
    
    y2 = y2 + rho2 * (Tu - z2);

    % === z3-update (enforce unit modulus) ===
    % tmp = u + y3/rho3;
    % z3 = tmp ./ abs(tmp);
    % y3 = y3 + rho3 * (u - z3);

    % === Dual updates ===
    

end

%% Final result
u_final = u;

z2 = T * u_final;
osu = ifft(cat(1, fft(u), zeros(3*N, 1)));
acf = conv(conj(flip(osu)), osu);
acfmag = abs(acf);
acfdb = 20*log10(acfmag / max(acfmag));
figure(1)
plot(acfdb, 'DisplayName', '优化后的自相关函数')

figure(2); plot(abs(T * u_final), 'DisplayName', '优化后的过采样包络');
% figure(2); plot(abs(osu), 'DisplayName', '优化后的过采样包络');

figure(3);
plot(abs(fft(T * u_final)), 'DisplayName', '优化后');

figure(4)
plot(acfmag, 'DisplayName', '优化后')






